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Abstract 



Neuronal activity arises from an interaction between ongoing firing generated sponta- 
neously by neural circuits and responses driven by external stimuli. Using mean-field 
analysis, we ask how a neural network that intrinsically generates chaotic patterns of ac- 
tivity can remain sensitive to extrinsic input. We find that inputs not only drive network 
responses, they also actively suppress ongoing activity, ultimately leading to a phase tran- 
sition in which chaos is completely eliminated. The critical input intensity at the phase 
transition is a non-monotonic function of stimulus frequency, revealing a "resonant" fre- 
quency at which the input is most effective at suppressing chaos even though the power 
spectrum of the spontaneous activity peaks at zero and falls exponentially. A prediction of 
our analysis is that the variance of neural responses should be most strongly suppressed at 
frequencies matching the range over which many sensory systems operate. 



Circuits of the central nervous system exhibit temporally irregular ongoing activity that is 
not directly related to sensory or behavioral events. The fact that this spontaneous activity 
is not suppressed by averaging over the large number of synaptic inputs to each neuron [1] 
suggests that chaotic network dynamics may represent a substantial local source of fluctuat- 
ing activity in cortical and subcortical circuits. Previous modeling studies have shown that 
nonlinear random network models with strong recurrent excitatory and inhibitory connec- 
tions generically exhibit chaotic dynamics [2, 3,4]. In this work, we ask how intrinsically 
generated fluctuating activity affects neuronal responses to external stimuli. The nonlin- 
ear effects of oscillatory drive, including frequency dependence and phase locking, have 
been well explored in low-dimensional chaotic dynamical systems (see e.g. [5, 6, 7, 8, 9]). 
However, relatively few studies have explored entrainment of extended, high-dimensional 
spatiotemporal chaotic systems by external forcing (see e.g. [10, 11, 12, 13, 14]). Here, we 
explore the locking of large chaotic neuronal networks to external stimuli and study how it 
depends on stimulus amplitude and frequency. 

We study phenomenological firing-rate network models representing neurons in a localized 
circuit that are coupled by relatively strong excitatory and inhibitory connections randomly 
distributed in the network. Specifically, we consider a network of N interconnected neu- 
rons, each described by an activation variable Xi for i — 1, 2, . . . N, satisfying 



with 4>{xi), which is a saturating monotonic function of the total synaptic input x^ repre- 
senting a normalized firing rate relative to a fixed background rate, r . Here we choose 



so that the normalized firing rate varies from zero to 2. For r = 1, we recover the often-used 
tanh function, but we use a smaller value of r = 0. 1, which is more biologically reasonable 
[15]. The time variable in Eq. 1 is defined in units of the single-neuron time constant, 
r r = 10 ms. Each element of the network connectivity matrix J is chosen randomly and 
independently [16] from a Gaussian distribution with zero mean and variance g 2 /N, where 
the gain g acts as the control parameter of the network. The external input term is set to 
Hi — I cos(ut + 9i), with the phase 9i chosen randomly and independently for each neuron 
from a uniform distribution between and 2n. This corresponds to situations in which 
the oscillatory input does not introduce global temporal phase coherence, which occurs, 
for example, for a population of neurons with a broad range of preferred spatiotemporal 
phases. 

To characterize the activity of the network, we make extensive use of the autocorrelation 
function of each neuronal rate averaged across all the units of the network, 
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Figure 1: Activity of typical network units (left column), average autocorrelation function (middle 
column) and log-power spectrum (right column) for a network with N = 1000, 5 = 1.5 and vq = 0.2. 
a) With no input (I = 0), network activity is chaotic, b) In the presence of a weak input (/ = 
0.04, / = uj/2tt = 4 Hz), an oscillatory response is superposed on chaotic fluctuations, c) For a 
stronger input (/ = 0.2, / = 4 Hz), the network response is periodic, d, e, f) Average autocorrelation 
function and g, h, i) Logarithm of the power versus frequency for the network states corresponding 
to panels a, b, and c. 

where the angle brackets denote a time average. C(0) is related to the total variance in the 
fluctuations of the firing rates of the network units, whereas C(r) for non-zero r provides 
information about the temporal structure of network activity. 

Previous work [2] has shown that, in the limit N — > oo with no input (/ = 0), this model 
displays only two types of activity: a trivial fixed point with all x = when g < 1 and chaos 
when g > 1 . The spontaneously chaotic state is characterized by highly irregular firing rates 
(Fig. la), a decaying average autocorrelation function (Fig. Id), and a continuous power 
spectrum (Fig. lg). Note that the fluctuations in Fig. la are considerably slower than the 
10 ms time constant of the model. The associated average autocorrelation function decays 
to zero as r increases (Fig. Id) implying that the temporal fluctuations of the spontaneous 
activity are uncorrelated over large time intervals, a characteristic of the chaotic state. The 
power spectrum decays from a peak at zero (Fig. lg) and, although it is broad, the power 
at high frequency is exponentially suppressed. Strong suppression of high-frequency fluc- 
tuations is another characteristic of the chaotic state in these networks. By comparison, the 
power spectrum of a non-chaotic network responding to a white-noise input falls off only 
as a power law at high frequencies. 

When this network is driven with a relatively weak sinusoidal input (Fig. lb, e & h), the 
single-neuron response consists of periodic activity induced by the input superposed on a 



3 



chaotic background (Fig. lb). The average autocorrelation function for the network driven 
by weak periodic input consequently reveals a mixture of periodic and chaotic activity (Fig. 
le). Periodic oscillations at the input frequency appear at large values of r, but the variance 
given by C(0) is larger than the height of the peaks in these oscillations. This indicates that 
the total firing-rate variance is not completely accounted for by the oscillatory response 
of the network to the external drive, the additional variance arising from residual chaotic 
fluctuations. Similarly, the power spectrum shows a continuous component generated by 
the residual chaos, a prominent peak at the frequency of the input, and peaks at harmonics 
of the input frequency arising from network nonlinearities (Fig. lh). 

When the amplitude of the input is increased sufficiently, the single-neuron firing rates 
oscillate at the input frequency in a perfectly periodic manner (Fig. lc), yielding a periodic 
autocorrelation function (Fig. If). C(0) now matches the height of the peaks in each of 
its subsequent oscillations, meaning that the periodic component in C accounts for the 
entire response variance quantified by C(0). All of the network power is focused at the 
frequency of the input and its harmonics, also indicating a periodic response free of chaotic 
interference (Fig. li). 

To explore these results analytically and more systematically, we developed dynamic mean- 
field equations appropriate for large N. The mean-field theory is based on the observation 
that the total recurrent synaptic input onto each network neuron can be approximated as 
Gaussian noise [2]. The temporal correlation of this noise is calculated self-consistently 
from the average autocorrelation function of the network. We begin by writing Xi — x^+x}, 
where x° is the steady-state solution to dx\jdt = —x® + I cos(cut + 9i) and x\ satisfies 
dxj/dt = —x\ + J2j Jij4>( x ] + x< j)- This implies that x®(t) = hcos(cut + Qj), where h = 
I/Vl + <j0 2 and we have incorporated a frequency-dependent phase shift into the factor 
Qi. Mean-field theory replaces the network interaction term in the equation for x\ by a 
Gaussian random variable 77, so that dxj/dt = —x\ + rji. Averages over time and network 
units (denoted by square brackets) as in Eq. 3, are implemented by averaging over J, 9 and 
77, an approximation valid for large N. 

Self-consistence is obtained in the mean-field theory by requiring that the first two moments 
of 77 match the moments of the network interaction that it represents. Thus, we set [r)i(t)] = 
J2j Jij4>{ x j{t)) — 0> because [J^] = 0. Similarly, using the identity [JaJjk] = g 2 8ij5 k i/N, 
we find that 
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Next, defining A(r) = [x\(t)x\(t + r)] and recalling that dx\/dt = —x\ + rji, it follows 
that 



d 2 A(r) 
dr 2 



= A(r) - g 2 C(r) . (5) 
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The final step in the derivation of the mean-field equations is to note that because x 1 (t) and 
x 1 (t+r) are driven by Gaussian noise, they are Gaussian random variables with moments 
[x 1 ^)] = [x l {t + t)\ =0 , [x 1 ^ 1 ^)] = [x 1 (t+T)x 1 (t+r)]=A(0), and [x^t+rfx 1 ^)] = 
A(r). To realize these constraints, we introduce three Gaussian random variables with zero 
mean and unit variance, zi for % — 1, 2, 3, and write 

x\t) = y/A(0)-\A(r)\ Zl + sgn(A(r))^\Mrj\z 3 



and 



x\t + r) = ^A(0)-\A(t)\z 2 + A(t)J\A(t)\z 3 . 



C can then be computed by writing x = x° + x 1 and integrating over these Gaussian vari- 
ables, 

C(r) = /^/^/^^(VAW-IAMI^ 
+sgn(A(r))^|A(r)|2;3 + /icos (0) 
x 1 1)^(^(0)- 1 A(r) |^ 2 + v/|A(r)|z 3 

+h cos (wr+0) ) , (6) 

where Dzj = dzjexp(— zf/2)/y/2ir for i = 1,2,3 and 6 = 9 + cot. Eq. 6 determines 
C(r) as a nonlinear function of A(r). Substituting this expression into Eq. 5 provides a 
nonlinear differential equation for A(r), with g, /i, w and A(0) as parameters. 

Eq. 5 has the form of the equation of motion for a classical particle of unit mass and position 
A(r) moving under the influence of a force that depends on C. This force is a function 
of the current position of the particle, A(r) (as well as on its initial position A(0)), and it 
contains terms representing external forcing that are periodic in r with period 2n/u. For 
weak inputs and g greater than but close to 1, Eq. 5 reduces to an undamped forced Duffing 
oscillator, although we do not restrict our analysis to this limit. 

The analogous mechanics problem has to be solved with the initial condition A(0) = 0, 
which imposes a smoothness constraint on the correlation function. The initial value A(0) 
is fixed by requiring that A(0) > A(r). We solved Eq. 5 numerically using iterative meth- 
ods to determine A(0), and found two types of solutions. The first is a solution in which 
A(t) is a periodic function of r with frequency to, as in Fig. If. This solution, which rep- 
resents a network state that is fully entrained by the oscillatory input, exists for all values 
of /, ijj and g. The second solution is characterized by A(r) that decays for small r and 
oscillates for large r, so that A(0) is larger than the peaks in the large-r oscillations, as 
in Fig. le. This solution, which corresponds to a non-periodic state only partially locked 
to the oscillatory drive, only exists for / smaller than a critical value that depends on u 
and g. A linear perturbation analysis of the mean field theory shows that this non-periodic 
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Figure 2: Phase transition curves showing the critical input amplitude that divides regions of pe- 
riodic and chaotic activity as a function of input frequency, a) Transition curves for tq = 0.2 and 
g = 1.5 (dashed) or g = 1.8 (solid). The stars indicate parameter values used in Figs, lb, e, h and 

1 c, f, i. The inset traces show representative single-unit firing rates for the regions indicated, b) 
A comparison of the transition curve computed by mean-field theory (open circles and line) and by 
simulating a network (filled circles) for tq = 1, g = 2 and, for the simulation, N = 10,000. 

solution is stable throughout the regime where it exists. The periodic solution is unstable 
in this regime and is stable outside it. The mean-field analysis also shows that the non- 
periodic solution corresponds to a state with "exponential" sensitivity to initial conditions 
(a positive Lyapunov exponent) [2], i.e., a chaotic state. 

The resulting phase diagram marks the transition between the periodic and non-periodic 
states (Fig. 2). Surprisingly, the transition curves are non-monotonic functions of frequency 
and reveal a "resonant" frequency at which it is easiest to entrain the chaotic network with 
a periodic input (even through there is no peak in the power spectrum of the chaotic activity 
at this frequency). This frequency is roughly twice the inverse time constant of the chaotic 
fluctuations in the spontaneous state and for g not too much greater than 1, the correspond- 
ing period can be an order of magnitude longer than the single-neuron time constant. Figs. 

2 & 3b indicate that internally generated fluctuations are most easily suppressed by stimuli 
oscillating in the few Hz range. 

The phase transition curve shifts upward and to the right as g increases (Fig. 2a & b), 
indicating a higher resonant frequency as well as a larger critical input amplitude. This 
occurs because the chaotic activity for larger g has a higher amplitude, making it more 
difficult to suppress, and a smaller inverse correlation time, leading to a higher resonance 
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Figure 3 : Signal and noise amplitudes as a function of input amplitude and frequency, a) Definition 
of the signal and noise amplitudes, a 2 sc and <7 c haos respectively, in terms of the mean-subtracted 
correlation function, b) Signal and noise amplitudes for / = 20 Hz, g = 1.5 and ro = 0.2 as a 
function of input amplitude. The transition from chaotic to non-chaotic regimes occurs at 7 = 0.44. 
c) Same as panel b, but with fixed input amplitude (I = 0.2) and varying input frequency. In the 
region between 3 and 7 Hz, responses of the network are free from chaotic noise. In b and c, open 
circles denote the signal amplitude and filled circles the noise amplitude. 

frequency. The location of the phase transition computed by mean-field theory is in good 
agreement with simulation results for large networks (Fig. 2b). 

To study the implications of the phase transition further, we divide network responses into 
signal and noise components by separating the full response variance into two terms, a 2 sc 
and of haos . For this purpose, we subtract the square of the average value of <p from C(r) and 
consider the mean- subtracted correlation function, C(r) — [(f)] 2 . The signal amplitude, cr osc , 
is the square root of the amplitude of the oscillatory part of this correlation function for 
large r (Fig. 3a). The noise amplitude, <r c haos, is the square root of the difference between 
the value of the mean- subtracted correlation function at r = and the peak of its oscillations 
(Fig. 3a). In the frequency domain, a 2 sc measures the total power in the network activity at 
the input frequency and its harmonics, whereas of haos measures the residual power. 

The signal amplitude increases linearly with the strength of the input (J) over the range 
considered in Fig. 3b. The noise amplitude has a more complex nonlinear dependence, 
reflecting the presence of the phase transition in Fig. 2 and duplicating the effect seen in 
Fig. 1, in which a sufficiently strong input completely suppresses the chaotic component 
of the response. An interesting feature to note is that there is no clear signature of this 
chaotic-to-periodic transition in the signal amplitude. When plotted as a function of input 
frequency for fixed /, the signal amplitude shows relatively weak frequency dependence 
below about 4 Hz and then rolls off at higher frequencies (Fig. 3c). This is a result of 
the low-pass filtering property of the network. The noise amplitude has a more interesting 
dependence. Between and 3 Hz, the noise amplitude drops steeply and vanishes for 
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frequencies between 3 and 7 Hz, rising again above 7 Hz. This double transition is a 
consequence of the non-monotonicity of the phase transition curves in Fig. 2. As in Fig. 
3b, there is no apparent indication of these transitions in the signal amplitude. 

It has previously been noted that chaotic activity in neuronal networks can be suppressed by 
either white-noise [13] or constant [14] input in discrete-time models. However, discrete- 
time versions fail to capture the rich dynamics of the chaotic fluctuations and their effect 
on responses to time-dependent inputs. Suppression of spatiotemporal chaos by periodic 
forcing has also been reported [10, 11, 12], mostly through numerical simulations. In some 
of these simulations, an optimal frequency for complete locking similar to Fig. 2 has been 
observed [10]. Our results show that such a resonance effect occurs even when the power 
spectrum of the unforced chaotic fluctuations falls monotonically from zero frequency (Fig. 
1). The networks we considered only describe the effects of fluctuations induced by local 
interactions, whereas additional sources of variability carried by long-range connections or 
by local sources of stochasticity are present in real neurons. Therefore, we predict that an 
experimental plot of response variability versus stimulus frequency will follow a non-zero 
U-shaped curve with a minimum in the several Hz range, rather than falling to zero as in 
Fig. 3c. 

Variability in cortical responses is sometimes described by adding stochastic noise linearly 
to a deterministic response [17, 18]. Our results indicate that the interaction between intrin- 
sically generated "noise" and responses to external drive is highly nonlinear. Near the onset 
of chaos, complete noise suppression can be achieved with relatively low amplitude inputs, 
weaker, for example, than the strength of the internal feedback. Thus, suppression of spon- 
taneously generated "noise" in neural networks does not require stimuli so strong that they 
simple overwhelm fluctuations through saturation. A number of experiments indicate that 
stimuli as well as attention can suppress firing-rate variability [19, 20, 21, 22, 23](but see 
[24]). Although other mechanisms for nonlinear suppression of neuronal variability have 
been proposed [25, 26, 27, 28, 29, 30], our analysis indicates that such suppression is a gen- 
eral property of the interaction between internally generated dynamics and external drive 
in a nonlinear network. 

Spontaneous fluctuations in neural activity occur across a wide range of timescales, with 
increasing variability over long time intervals [31] and increasing power at low frequencies, 
although resonances may appear [24, 32]. In this work we have focused on firing-rate fluc- 
tuations using smooth rate-based dynamics, not spiking dynamics. Spiking neuron models 
with strong 'balanced' interactions can exhibit chaotic firing patterns [2, 3], but the fluctu- 
ations they produce have relatively flat power spectra associated with variability in short 
interspike intervals. It will be interesting to study stimulus effects in spiking network mod- 
els that exhibit slow irregular modulations of firing rates. 

In our model, weak correlations (of the order of 1/ y/N) in activity fluctuations exist be- 
tween all pairs of neurons. These correlations are distributed evenly between negative and 
positive values across the population. Slow spontaneous rate fluctuations in the cortex are 
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often associated with long-range spatial correlations, especially in anesthetized animals 
[33, 34]. As in our model, the observed spatial correlations are weaker than the firing rate 
autocorrelations. In some cases, both negative and positive rate fluctuations are also ob- 
served, such that the mean value of the pairwise correlations across a populations is much 
smaller than the width of the distribution of correlations [35, 36, 37]. However, the extent 
of the contribution of local network dynamics to the observed low frequency correlations 
is unclear [22,34]. 

Neuronal selectivity to stimulus features is typically studied by determining how the mean 
response across experimental trials depends on various stimulus parameters. The presence 
of nonlinear interactions between stimulus -evoked and spontaneous fluctuating activity in- 
dicates that response components that are not locked to the temporal modulation of the 
stimulus may also be sensitive to stimulus parameters. In general, our results suggest that 
experiments studying the stimulus-dependence of the noise component of neural responses 
could provide important insights into the nature and origin of activity fluctuations in neu- 
ronal circuits, as well as their role in neuronal information processing. 
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